/*
 * AqueousPhase.h
 *
 *  Created on: 21 Jul 2011
 *      Author: al2010
 */

#ifndef AQUEOUSPHASE_H_
#define AQUEOUSPHASE_H_

// Eigen includes
#include <Eigen/Dense>
using namespace Eigen;

// GeoReact includes
#include "AqueousActivity.h"
#include "Phase.h"
#include "Utils.h"

class AqueousPhase : public Phase
{
public:
	enum IonicActivityModel { Davies, ExtendedDebyeHuckel, Hkf, Pitzer };
	
public:
	AqueousPhase();
	
	AqueousPhase(const StringVector& speciesNames);
	
	virtual ~AqueousPhase();
	
	/// --- Setter methods --- ///
	void
	SetSpecies(const StringVector& species);
	
	void
	SetIonicActivityModel(IonicActivityModel model);
	
	void
	SetSetschenowConstant(String speciesName, const DoubleVector& bi, const DoubleVector& Ti);
	
	/// --- Getter methods --- ///
	const StringVector
	GetIonicSpecies() const;
	
	const DoubleVector
	GetIonicCharges() const;
	
	const Index
	GetWaterIndex() const;

	/// --- Computation methods --- ///
	const double
	EffectiveIonicStrength(const DoubleVector& n) const;
	
	const double
	StoichiometricIonicStrength(const DoubleVector& n) const;
	
	const double
	TotalMolalitySolutes(const DoubleVector& n) const;
	
	const DoubleVector
	Molalities(const DoubleVector& n) const;
	
	const DoubleVector
	IonicMolalities(const DoubleVector& n) const;
	
	const DoubleVector
	IonicStoichiometricMolalities(const DoubleVector& n) const;
	
	/// --- Phase virtual methods --- ///
	virtual	const String
	GetName() const;
	
	virtual const StringVector
	GetSpecies() const;
	
	virtual const DoubleVector 
	ActivityCoefficients(double T, double, const DoubleVector& n) const;
	
	virtual const DoubleVector
	Activities(double T, double, const DoubleVector& n) const;
	
private:
	/// The index of the water species.
	Index iwater;
	
	/// The vector of electrical charges of the ionic species.
	DoubleVector zi;
	
	/// The indexes of the ions, neutral species and electrolytes.
	IndexVector ions, neutrals, electrolytes;
	
	/// The stoichiometric matrix between the electrolytes and ions.
	MatrixXd vei;
	
	/// The function that computes the activity coefficients of the ionic species.
	IonicActivityFunction ionicActivityFunc;
	
	/// The function that computes the activity coefficients of the neutral species.
	vector<NeutralActivityFunction> neutralActivityFuncs;
	
private:
	void 
	InitElectrolytes();
	
	void 
	InitDefaultActivityModels();
	
	const DoubleVector
	ActivityCoefficientsIonicSpecies(double T, double P, const DoubleVector& n) const;
	
	const DoubleVector
	ActivityCoefficientsNeutralSpecies(double T, double P, const DoubleVector& n) const;
};

#endif /* AQUEOUSPHASE_H_ */
